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Abstract 

Pencils of Hankel matrices whose elements have a joint Gaussian distribution 
with nonzero mean and not identical covariance are considered. An approx- 
imation to the distribution of the squared modulus of their determinant is 
computed which allows to get a closed form approximation of the condensed 
density of the generalized eigenvalues of the pencils. Implications of this re- 
sult for solving several moments problems are discussed and some numerical 
examples are provided. 
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Let us define the random Hankel matrices 



do cti • ■ • (^p—l 
Ul 0,2 ■■ ■ dp 
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(32 0-3 
Clp flp+l 



Op 



(1) 



0, 1, 2, . . . , n — 1, efc is a complex 



ap_i Op ... a„_2 

where n = 2j9, 0^ = 5^ + e/c, 
Gaussian, zero mean, white noise, with variance cr^ and Sk G ■ 

Let us consider the generahzed eigenvalues {^j, j = 1, . . . ,p} of (f/i, Uq) 
i.e. the roots of the polynomial P{z) = det[Ui — zUq)] and the associated 
condensed density h{z), introduced in [8], which is the expected value of the 
(random) normalized counting measure on the zeros of P{z) i.e. 

1 



h(z) = -E 
P 



.i=i 



or, equivalently, for all Borel sets A cCT 

p 



J h{z)dz = -Y^^ProbiQ E A). 



It can be proved that (see e.g. [Ij) h{z) = -^Au{z) where A denotes 
the Laplacian operator with respect to x, y if z = x + iy and u{z) = 

{log{\P{z)\'^)} is the corresponding logarithmic potential. 

The condensed density h{z) plays an important role for solving moment 
problems such as the trigonometric, the complex, the Hausdorff ones. It was 
shown in [9l Uni El El d HJ |2l [3] that all these problems can be reduced to the 
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complex exponentials approximation problem (CEAP), which can be stated 
as follows. Let us consider a uniformly sampled signal made up of a linear 
combination of complex exponentials 

s. = XIc4. (2) 

where Cj,^j G (T. Let us assume to know an even number n > 2p of noisy 
samples 

ak = Sk + ek, k = 0,l,2,...,n-l 

where is a complex Gaussian, zero mean, white noise, with finite known 
variance o"^. We want to estimate p,Cj,^j, j = l,...,p, which is a well 
known ill-posed inverse problem. We notice that, in the noiseless case and 
when n = 2p, the parameters C,j are the generalized eigenvalues of the pencil 
(f/i, Uq) where now Uq and Ui are built as in ([T]) but starting from {sk}- 

From its definition it is evident that the condensed density provides infor- 
mation about the location in the complex plane of the generalized eigenvalues 
of {Ui, Uq) whose estimation is the most difficult part of CEAP. Unfortunately 
its computation is very difficult in general. In [7] a method to solve CEAP 
was proposed based on an approximation of the condensed density. An ex- 
plicit expression of h{z) proposed by Hammersley [8] when the coefficients of 
P{z) are jointly Gaussian distributed was used. The second order statistics 
of these coefficients in the CEAP case were estimated by computing many 
Pade' approximants of different orders to the Z-transform of the data {flfc}- 
This last step was essential to realize the averaging that appears in the defini- 
tion of h{z), which is the key feature to make the condensed density a useful 
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tool for applications. In fact in the noiseless case h{z) is a sum of Dirac 5 
distributions centered on the generalized eigenvalues while, when the signal 
is absent {s^ = VA;), it was proved in [Ij that if z = re*^, the marginal 
condensed density h^^'\r) w.r. to r of the generalized eigenvalues is asymp- 
totically in n a Dirac 6 supported on the unit circle Vcr^. Moreover for finite 
n the the marginal condensed density w.r. to 6 is uniformly distributed on 
[— 7r,7r]. Therefore if the signal-to-noise ratio (e.g. SNR = ^ min/j=i p |c/i|) 
is large enough h{z) has local maxima in a neighbor of each j = 1, . . . ,p 
and this fact can be exploited to get good estimates of ^j. However usually 
we have only one realization of the discrete process {a^}, hence we cannot 
estimate h{z) by averaging. We then look for an approximation of h{z) which 
can be well estimated by a single realization of {ofc}- The specific algebraic 
structure of CEAP will be taken into account and it will be shown that the 
noise contribution to h{z) can be smoothed out to some extent simply acting 
on a parameter of the approximant. 

The paper is organized as follows. In Section 1 some algebraic prelimi- 
naries are developed. In Section 2 the closed form approximation of h{z) is 
defined. In Section 3 we show how to get a smooth estimate of h{z) from the 
data by exploiting its closed form approximation. Finally in Section 4 some 
numerical examples are provided. 
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Let us consider the pencil F = Ui ~ zUq. The problem can be reduced to 
real (random) variables by using the following result (jSl Theorem 5.1]): 

Proposition 1 If F = Ui - zUo, z eCT and F = Vr + iVj, Vr, Vj G 
then 

\det{F)\^ = det{G) 



where G is the real isomorph of F, i.e. 



G = n{F) 



Vr -Vi 

Vi Vr 



Let us consider the QR decomposition of G where Q^Q = QQ^ = In 
where T denotes transposition, R is an upper triangular matrix and /„ is the 
identity matrix of order n. We then have 

\det[Ui - zUo)]\'^ = det{G) = \det{G)\ = \det{QR)\ = \det{R)\ = JJ \Rkk\. 

k=l,n 

We are therefore interested on the distribution of |-Ra;A:|, k = 1, . . . ,n. In order 
to perform the QR decomposition of the random matrix G = [g_^, . . . , g^] we 
make use of the Gram-Schmidt algorithm because it produces a triangular 
matrix R with positive elements. It is given by where Q = [q^, . . . ,q^]. 
For A; = 1 we get 
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ror K — i, . . . , /i 








if /c > 1 then 




^ik — 1^ y_^-i * - 


- 1 i- — 1 

- 1, . . . , /v 1 




1 -Rifc^j 


end 




Rkk = V^Zfe^^jt 




wife 

'lk~ 




end 





where, denoting by 3?, 3" the real and imaginary part of a complex number, 
= [3?(ao), . . . , 3?(ap), ^^(ao), • • • , ^(ap)]^ 



and 



^7 



Zr 



1 
-z 1 








1 



In fact by using the property of the real isomorph 

n{A)[a,bf = [^{Ac),'^{Ac)f, c = a + ib 
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if = [a^ji, tti/]"^ and is the k-th column of /„, we have 

9, = G[e,,Of = n{F)[e,,Or = mFe,),^{Fe,)r 
= [^[{Zn + iZj){a^ji + ia^j)], ^[{Zr + iZj){a^R + ia-^j)]f 

= n [{Zr + iZj)] [a^R, a^jf = Za^. 

We notice that is a random Gaussian 2(p + 1) vector with mean 



^2 



and covariance matrix S = ^hip+i)- Therefore Ru is distributed as the 
square root of a quadratic form in normal non-central variables. This is a 
well known distribution (see e.g. pT]). 
For k = 2 we have 



w 



919 



-2 -1 -2 £^^^-1 

and 



-R22 = Y ^2^2. 

Therefore R22 is no longer distributed as a quadratic form. However we 
notice that the vector W2 is a linear combination of the Gaussian vectors 

T 

and q with weights 1 and =7^. If this random weight were deterministic, 
R22 would be distributed as the square root of a quadratic form in normal 
non-central variables. We then consider the following approximation. We 
first notice that in the Gram-Schmidt algorithm, w^, k = 1, . . . ,n can be 
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computed as 



w 



1 = £i 



fc-i 



k-l 



Wv 



Ik-^Rikg^^ 9f.-^WikWi, Wik^Rik/Rii, k^2,..., 



i=l 



We then apply the Gram-Schmidt algorithm to the matrix 



S = E[G] 



E[Vr] E[Vi] 
-E[Vj] E[Vr] 



where 



n 



E[Vr] + iE[Vi] = E[F] =S,- zSo 

and So, Si are the Hankel matrices built from {sk} = {E[ak]} in the same 
way as Uq and Ui are built from {ak}- Let S = QR the QR decomposition 
of the deterministic matrix S, and consider the approximation of the vectors 
m given by: 



Wi = 



k-l 



Wi 
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WjjfcWj, Wik = Rik/Rii, k^2,...,n 



i=l 



It turns out that Wj^ are linear combinations of the Gaussian vectors g^,...,g^ 
with deterministic weights 1 and Wik and we will prove in the next section 
that 



Rkk ~ \/wlwk 
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which is distributed as the square root of a quadratic form in normal non- 
central variables. In fact we have 







In 
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W^,n 
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W2,n-In . 
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{WMQw (3) 



where 



W = 



1 

1^12 1 



... 
... 



eiR"^", g = {In6dZ)a, anda = 



Wi,n W2,n ... 1 

with ttfc e 1R^+^ and 

[3?(afc_i), . . . , K(afc+p_i), $5(afc_i), . . . , $5(afc+p_i)]^ 



an 



if /c = 1, . . . ,p 

[-Q^(afe_i), . . . , -Q^(ajk+p_i), K(ajk_i), . . . , ^{ak+p-i)Y if /c = p + 1, . . . , n 
hence, if — e,j^^ In and e^. is the k—ih. column of the identity matrix 7„, 
we have 



f{W-^ 60 4) (e,e^ 60 In) {W-^ 60 /n)£ = A/£^(W^-^e,e^W^-i (SO 4)^ = 



a^(/„ (SO Z^)(iy-^e,e^iy-i (SO (SO Z)a = A/a^(M^-^e,e^iy-i (SO ZTZ)a. 



But we notice that if a = [3?(ao)) ■ ■ ■ , ^(O'n-i), ^(^o)) ■ ■ ■ , ^{O'n-i)]'^ then a — 
Ba where 



fc-i 



Bfe = [C..,o,Vi,o,...,o]eiR(^+^)^ 
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' Bf . . . . -Bf 

B^ ^ 

Bf . . . . B^ Bf 

Summing up, if Ak = B^{W-^e^^elW-^ (g) Z^Z)B e iR2nx2n ^^^^^ 
shown that 

where a ~ N{s,^l2n) and s = [K(so), ■ ■ ■ , 3f^(sn-i), ^^(•so), ■ ■ ■ , Q'(sn-i)]^]- 
Moreover the following Proposition holds true 

Proposition 2 If Ak — DkPk is the spectral decomposition of A^ e 
R'^^x^^, then 

a^AkO, = alAkttk, 

where Dk = diag[Ak,0], Afc G R"^"", = [4, 0„xn]-Pfea, and Ak > 0. 
Moreover Uk ~ ^(/f^, ^--^n) w/iere /x^ = [/„, 0„xn]-PfcS- 

Proof. By construction we know that 

n 

a^AkU^wlwk = X^tL 

h=l 

and Ak > 0. Hence by the Sylvester's law of inertia n — rank{Ak). Therefore 
the diagonal matrix Dk of the eigenvalues of Ak can be partitioned as Dk — 
diag[Ak,0\ where > 0. Letting x = PkO, and = [a^,?/"^], we have 

X = alAk&k. 
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Bl 



T 



jp^n{n+2)x2n 



a' Aka = a^P^DkPka = a' Pk 



Tr>T 



Ak 




DkPka = X 



Ak 
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The last part of the thesis follows by noticing that //^ = E[a^] — [7„, Onxn\PkE[a\ 
and 

E[{a,-E[a,]){a,-E[a,]) = [In,On><n]PkE[{a-E[a]){a-E[a]f]P^[In,On><nf = 

□ 

2 Closed form approximation of h{z) 

In order to build a closed form approximation of the distribution of u{z) we 
notice that we can rewrite u{z) as 

u{z) = Ie {log{\P{z)\')} = ^E {log{\det[F]n = 

lE{\og{\det[Gf)]^lY.E[\og{Rl,)\ 

lb lb 

k=l 

Hence we will consider in the following the distribution of i?^^ and its ap- 
proximation by the quadratic form in normal variables O'^^jta. 

Lemma 1 For k — 2, . . . ,n; i — 1, . . . ,k — 1 let be Wik — Rik/Ru where 
G = QR is the QR decomposition of the random matrix G and Wik = Rik/Ru 
where S = QR is the QR decomposition of the matrix S = E[G]. Then, for 

(7^0, 

E[{Wik-mkf]^o{a). 
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Proof. Let us consider the rational function fik{a) — Wik and its Taylor 
expansion around s — E[q\: 



fik{a) = fik{s) 
which can be rewritten as 



and, squaring and taking expectations, remembering that a ~ -/V(s, ^/2n) 
and applying Holder's inequality we get 

2" 



E[{Wik-Wikf] = E 



\ h J 



+ ... 



hk 



hk 



21 \ 1/2 



h^^k 
h 

h^k h 

The thesis follows because the dropped terms are functions of {cr^^, h > 2}. 
□ 

Theorem 1 Let x and x be the vectors obtained by stacking respectively the 
elements Wik, o-nd Wik, for {/c = 2, . . . , n; i = 1, . . . , A; — 1}. Let be qk{x) — 
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Rli^ where G = QR is the QR decomposition of G, and qk{x_) = a'' A^a where 
Ak = B^{W-^e^elW-^ (g) Z^Z)B. Then, for a ^ 0, 

E[\qk{x) - qk{x)\] = o{v^). 

Proof. Let us consider the rational function qk{x) and its Taylor expansion 
around x: 

qk{x) = qk{x) + ^ \a:=x {Xh - Xh) + . . . . 

^-^ OXu ~ ~ 

h —'^ 

But, by definition, quix) = qk{x), hence, squaring and taking expectations, 
we have 

E[\qk{x) - qk{x)\] = E[\J2Ph{xh-Xh)\] + ... 

h 

< E[\Ph{xh — Xh)\] + . . . by Minkowsky's inequality 

h 
h 

by Holder's inequality. The thesis follows from Lemma [1} □ 

Theorem [l] allows us to assume that {Rlk, k = 1, . . . ,n} are approxi- 
mately distributed as a quadratic form in normal variables provided that the 
noise is small enough. A large literature exists about distribution of quadratic 
forms in normal variables (see e.g. [H] and references therein). Several series 
expansions are considered for the density of the quadratic form. However we 
will make use of the approximation proposed in [13] which makes it possible 
to get a closed form expression of E[\og{Rl,^)]. 
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Theorem 2 Let be aj^ ~ ^(l^k^ T-^") '^^^ ~ ^fc^fe^fc ^ quadratic 

form with rank{Ak) = n. Then there exist positive constants ak,(3k,1k such 
that the first three moments of jSkiXaf.V'' rnatch those of q/^. 

Proof. Let us denote by Uhk, h = 1,2,3 the first three moments w.r. to the 
origin of g^, which exist and whose explicit expression can be found in pT| 
Theorem 3. 2b. 2]. The corresponding moments of (3k{Xa 



2^^ r(7fc+afc/2) 



r(afc/2) 

(2/?fc 



(nn N2 r(27fc+afc/2) 

l^Pfc; r(afc/2) 



3 r(37fc+afc/2) 
r(afc/2) 



Then is obtained by solving 

ra + «fc/2) 



i/ifc, i.e. 



z/ifcr(dfe/2) 



r(d,/2) • ■ "'^ 2r(7, + dfc/2) 

where a?; and 'jk are obtained by solving the system of equations 



r(afc/2)r(27fc+afc/2) 
[r(7fc+afe/2)]2 

[r(afc/2)]2r(37fc+«fc/2) 
[r(7fe+afc/2)]3 



which admits an unique positive solution. In fact the mapping 

r(a/2)r(27+a/2) 



/(a, 7) = 

is continuous in M'^ x iR+ and 



[r(7+a/2)]2 

[r(«/2)]2r(37+a/2) 

[r(7+a/2)]3 



2/i(«, 7) [^(a/2 + 27) - + 7)] 

3/2(a, 7) + 87) - ^(«/2 + 7)] 



> 
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because /(a, 7) > and = is monotonic increasing for x E IR^. 

Therefore fi{a, 7) and /2(tt, 7) are monotonic increasing for each fixed a > 
as functions of 7. Moreover we notice that f(a,0) = l,Va; and di = ^ > 
1; (^2 = ^ > 1 by Holder inequahty, taking into account the positivity of 
the r.v. gfc. Therefore /(a, 7) can assume every value in [1, 00). □ 
We can now prove the main theorem 

Theorem 3 For cr ^ 0, and approximating the distribution of Rlf.{z),k = 
1, . . . ,n by the density Pkiz){x1.^(^^-^)'^''^''\ where Q{z)R{z) is the QR decom- 
position of G{z), which is the real isomorph of Ui — zUq, we have 

u{z) = ^E{log{\det[Ui - zUo]]"^)} ^ u{z) 

where 

1 " 

^(^) = - 5Z[l°g(2)7fc(^) + log[/3,(z)] + ^kiz)nMz)/2] 
k=i 

and 

h{z) ^ ^Au{z). (4) 



Proof. By Proposition [l] the computation of \det[Ui — zUo]\'^ is reduced to 
that of Rli^, k = 1,. . . ,n, whose density, for each k, can be approximated 
for (7 ^ by that of a quadratic form in normal, noncentral, variables by 
Theorem [Tj By Proposition [2] this quadratic form is equivalent to a quadratic 
form with diagonal positive matrix whose density can be approximated, up to 
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the third moment, by the density PkixiJ^'' for suitable constants (yk,Pk,1k 
by Theorem [2] But, if g is a r.v. with density (3{Xa)^, we have 

1 r°° 

E[log{q)] = ^^i^^^^i^^ \ogW)q'^"-'e-^'^dq = log(2)7+log(/3)+7Vl>(a/2). 
□ 



3 Smooth estimate of the condensed density 

We want to show now that we can exploit the closed form expression of the 
condensed density to smooth out the noise contribution to h{z). This allows 
us to get a good estimate of p and ^j,j = 1, . . . ,p which is the nonlinear most 
difficult part of CEAP, from a single realization of the measured process {flfc}- 
We first notice that if in Theorem [2] we look for an approximant which 
fits the first two moments of then the solution is 



Ik = 1, Pk = ^ , ttfc 



^2 k - T^lk . _ '^^Ik 

o ' 2 

•^I^lfc l^2k — l^ik 



and ElRlf^] = vik = &kPk- Hence, dropping the terms independent of z and 
using the approximation ~ log(x) — - (see[T2]) we get 



1 " 

^(^) ^ — $^A(log[/3fc(;^)] + M/K(^)/2]) 



4?™ 

k=l 



1 " / 

-5^A log[/?fc(z)]+logK(z)/2] 



4™ V Mz) 



-X:AflogK(^)/2]-^y (5) 
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Now we notice that 



i^2k - 1^1, = ^ and = ./^ (6) 

Moreover from pTi Theorem 3. 2b. 2] the first two moments w.r. to the origin 



of are 



^2. = ^ELiA.W^ + 2aWu + 



'Ik 



where Xk{h) are independent of cr^ because they are the (positive) eigenvalues 
of Ak which is independent of cx^. Therefore when s 7^ if 

^ n n 



/i=i 



/i=i 



Pk = — z = cr 



2z^ifc 

xpresj 

measured e.g. by p — 



a^l + 26' 
a^c + 2d 



2 (^^c + c/ 



Deriving these expressions respectively with respect to cr^ and to the SNR 



But 



- ^ we get 

dpk _ bd + aa\2d + ca^) 
~ {d + ca^y 

^ _ _ Hl^f^Hk _ 



> 0, 



d&k d{c + p){2ad + b{p — c)) 



dp 



{ad + 6p)2 



p-c 



n n / 

h=\ h=l ^ 



ay 2 



- 1 Xkih) > 



if > ^ V/i. Hence is an increasing function of cr^ and cxk is an increasing 
function of p if the SNR is not too small consistently with the expression (|6|. 
When a^ = 

Pk = 0, akPk = I'lk = f^l^kP,, lim ak = 00 
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and, by (|5| 

1 " 



k=l 



The last equality can be checked also by noticing that 



The idea is then to use the parameters j3k as smoothing parameters and 



as signal-related parameters. By fixing = \/k and taking cxk — — 



the variance of R\k{z) is controlled by [3 and h{z) can be estimated by 

where A is the discrete Laplacian and vik{z) is an estimate of vik{z). In the 
following we use, as an estimate of I'lkiz), the value Rl^iz) obtained by the 
QR decomposition of the data matrix G{z) defined in Proposition 1. 

From a qualitative point of view increasing (3 has the effect to make larger 
the support of all modes of h{z) and to lower their value because h{z) is a 
probability density. Hence the noise-related modes are likely to be smoothed 
out by a sufficiently large p. However a value of f3 too large can result in a 
low resolution spectral estimate. 

4 Numerical results 

In this section some experimental evidence of the claims made in the previous 
sections is given. 
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To appreciate the advantage of the closed form estimate h{z) with respect 
to an estimate of the condensed density obtained by MonteCarlo simulation 
an experiment was performed. N = 100 independent realizations of the 
r.v. ci'k\k = l,...,n, r = 1,...,N were generated from the complex 
exponentials model with p = 5 components given by 

C _ r -0.1-i27r0.3 -0.05-i27r0.28 -0.0001+j27r0.2 -0.0001+i27r0.21 -0.3- j27r0.351 
S_ I 1 1 1 ' J 

c= [6,3,1,1,20], n = 74. 

We notice that the frequencies of the 3'''^ and 4*^^ components are closer than 
the Nyquist frequency (0.21 — 0.20 = 0.01 < 1/n = 0.0135). Hence a super- 
resolution problem is involved in this case. Two values of the noise s.d. a 
were used 

(J = 0.2, 0.8. 

An estimate of h{z) was computed on a square lattice of dimension m = 100 
by 

N n / 
r=l k=l \ 

where R^'^\z) is obtained by the QR factorization of the matrix G{z) built 
from sample r-th. In the top part of figjl] the estimate of h(z) obtained by 
Monte Carlo simulation is plotted. In the bottom part the smoothed esti- 
mates h{z) for a = 0.2 and (3 = 5n based on a single realization was plotted. 
In figj2]the results obtained with a = 0.8 and (3 = 5n are shown. We notice 
that by the proposed method we get an improved qualitative information 
with respect to that obtained by replicated measures. This is an important 



K 



kk y 



2(t2/5 
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feature for applications where usually only one data set is measured. We also 
notice that when a = 0.2 the probability to find a root of P{z) in a neighbor 
of is larger than the probability to find it elsewhere. This is no longer 
true when a — 0.8 even if the signal-related complex exponentials are well 
separated. In the following we will say that the complex exponential model 
is identifiable if this last case occurs and it is strongly identifiable if the first 
case occurs. Therefore if the model is identifiable the signal-related complex 
exponentials are well separated but the relative importance of some of them 
- measured by the value of the local maxima of h{z) - is not larger than the 
relative importance of some noise-related complex exponentials. Therefore 
in this case we need some a-priori information about the location of the 
in order to separate signal-related components from the noise-related ones. 

We want now to show by means of a small simulation study the quality of 
the estimates of the parameters p,^ and c which can be obtained from h{z). 
To this aim the following estimation procedure was used: 

• the local maxima of h{z) are computed and sorted in decreasing mag- 
nitude 

• a clustering method is used to group the local maxima into two groups. 
If the model is strongly identifiable the signal-related maxima are larger 
than the noise-related ones, therefore a simple thresholding is enough 
to separate the two groups. A good threshold is the one that produces 
an estimate of Sk which best fits the data Ofe in L2 norm as the noise is 
assumed to be Gaussian 
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• the cardinality p of the class with largest average value is an estimate 
of p 

• the local maxima — 1, . . . , p of the class with largest average value 
are estimates of ^j, j = 1, . . . ,p. Of course if p p some are not 
estimated or viceversa some spurious complex exponentials are found 

• c is estimated by solving the linear least squares problem 

c = argminx\\Vx — d\\l 
where V e (T"^^ is the Vandermonde matrix based on ^j, j = 1, . . . ,p 

The bias, variance and mean squared error (MSE) of each parameter sep- 
arately were estimated. = 500 independent data sets a^^^ of length n 
were generated by using the model parameters given above and a = 0.2. 
For r = 1, . . . , the condensed density estimate h^'^\z) was computed on 
a square lattice of dimension m — 100. The estimation procedure is then 
applied to each of the h^^\z),r = 1, . . . ,N and the corresponding estimates 
ij~\cj'\j — 1, . . . jP*^''-' of the unknown parameters were obtained. If the 
estimate p^^^ was less than the true value p, the corresponding data set a^'^^ 
was discarded. 

In Table 1 the bias, variance and MSE of each parameter including p is 
reported. They were computed by choosing among the — 1, . . . ,p^^^ > 

p the one closest to each ^k,k — 1, . . . ,p and the corresponding Cj . If more 
than one is estimated by the same i^j"^ the r— th data set a^''^ was discarded. 
In the case considered all the data sets were accepted. 
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As a second example the reconstruction of a piecewise constant function 
from noisy Fourier coefficients is considered. The problem is stated as follows. 
Given a real interval [— tt, tt] and + 1 numbers — tt < li < I2 ■ ■ ■ < /jv+i < tt, 
let be the class of functions defined as 

N 

where 

Xj{t) = < 

I otherwise 

and the wj are real weights. The problem consists in reconstructing a function 
F{t) ^ T from a finite number of its noisy Fourier coefficients 

1 r 

cifc = - y F{t)e'^^dt + Cfc = Sfc + efc , A; = 0, . . . , n - 1 , 

where efc is a complex Gaussian, zero mean, white noise, with variance . 
We are looking for a solution which is not affected by Gibbs artifact and can 
cope, stably, with the noise. The basic observation is the following. The 
unperturbed moments are given by 

'^^\\[ ^(^y^'dt = wj'^^^ exp(iA,fc), 



where 



3=1 



Pj - ^ ' ^3 



2 ' ^ 2 
Then consider the Z-transform of the sequence {sk) 



/ 1 z- e}^^ \ 
(z) = > Wi { 3i H In -, — 

.7=1 ^ ^ 
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which converges if 1^1 > 1 and is defined by analytic continuation if \z\ < 1. 
We notice that s{z) has a branch point at C,j = e*'^ , j = 1, . . . , + 1 where 
Ij are the discontinuity points of F{t). It was proved in [9l [10] that the Cj 
are strong attractors of the poles of the Pade' approximants [m,n]f{z) to the 
noisy Z-transform 

oo 

f{z) = Y,akZ-'' 

k=0 

when m,n —* oo and m/n — 1. It is easy to show that the poles of [m, nlflz) 
are the generalized eigenvalues of the pencil (f/i,f/o) built from the data 
ak, k = 0, . . . ,n — 1 whose condensed density is h{z). Therefore, as shown in 
[HI [To] the local maxima of h{z) are concentrated along a set of arcs which 
ends in the branch points C,j and on a set of arcs close to the unit circle. As 
the branch points are strong attractors for the Pade' poles, the probability to 
find a pole in a neighbor of a branch point is larger than elsewhere, therefore 
it can be expected that the branch points correspond to the largest local 
maxima of h{z), as far as the SNR is sufficiently large. In order to compute 
estimates of it is sufficient to compute the arguments of the main local 
maxima of h{z). The Wj are then estimated by taking the median of the data 
in each interval The median is in fact robust with respect to errors 

affecting Ij. 

The method was applied to an example considered in [TU] where compar- 
isons with other methods were also reported. In the top left part of figjS] 
the original function F{t) is plotted. In the top right the noisy data with 
SNR = 7 are reported where the SNR is measured as the ratio of the stan- 
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dard deviations of {s^} and {cfc}. In the bottom parts the condensed density 
and the reconstructed function F(t) are plotted. Looking at the condensed 
density we notice that the model is strongly identifiable, therefore the esti- 
mation procedure outlined above was applied. In fig j4] the same quantities as 
above but with SNR = 1 are plotted. In this case the model is identifiable 
but not strongly therefore the clustering step does not work. The number of 
complex exponentials used to get the reconstruction plotted in fig|4]is p = 20 
and was found by trial and errors. 

We notice that when SNR = 7 we get an almost perfect reconstruction, 
better than that reported in [TU] . When SNR = 1 the reconstruction quality 
is worse as expected but still comparable with the one reported in [TO] . 
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P 


bias(p) 


s.d.(p) 


MSE{p) 




5 


0.0000 


0.0000 


0.0000 






bias(^-i) 


s.d.^-i 


MSE(ii) 


J = l 


-0.2796 - 0.8606i 


-0.0008 + O.OOOli 


0.0000 


0.0000 


J = 2 


-0.1782 - 0.9344i 


0.0036 - O.OOlOi 


0.0000 


0.0000 


J = 3 


0.3090 + 0.9510i 


0.0057 - 0.0064i 


0.0031 


0.0001 


J =4 


0.2487 + 0.9685i 


-0.0058 + 0.0110 


0.0019 


0.0002 


J = 5 


-0.4354 + 0.5993i 


-0.0047 + 0.0054i 


0.0108 


0.0002 






bias{cj) 


s.d.{cj) 


MSE{cj) 


J = l 


6.0000 


0.0440 


0.1238 


0.0173 


J = 2 


3.0000 


-0.0407 


0.0688 


0.0064 


J = 3 


1.0000 


O.OIll 


0.0736 


0.0071 


J =4 


1.0000 


-0.6767 


0.0808 


0.4644 


J = 5 


20.0000 


-0.1007 


0.2574 


0.0764 



Table 1: Statistics of the parameters p, ^j, j — 1, . . . ,p and Cj,j — 1, . . . ,p 




Figure 1: Top: Monte Carlo estimate of the condensed density when a — 0.2; 
bottom: estimate of the condensed density by the closed form approximation 
with (5 = 14.8. 
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Figure 2: Top: Monte Carlo estimate of the condensed density when a — 0.8; 
bottom: estimate of the condensed density by the closed form approximation 
with P = 23.7. 
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Figure 3: Top left: original function; top right: original function with Gaus- 
sian noise added (SNR=7). Bottom left: estimate of the condensed density 
by the closed form approximation; bottom right: reconstruction of the orig- 
inal function. 
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Figure 4: Top left: original function; top right: original function with Gaus- 
sian noise added (SNR=1). Bottom left: estimate of the condensed density 
by the closed form approximation; bottom right: reconstruction of the orig- 
inal function. 



